lambda_b = 0.9;
lambda_g = 0.85;
gamma = 0.2; 
r_hat = 0.8;
theta_b = 1;
theta_g = 1;
theta_bb = 1;             %theta_B+
theta_gg = r_hat;         %tehta_G+  
p = 1;
q = 0.1;    
U = 5;


r = r_hat:0.005:1;
n = length(r);

x0 = [0,0,0,0];
A = [];
b = [];
Aeq = [];
beq = [];
lb = [];
ub =[];
options = optimoptions('fmincon','Display','iter');

X_sol = ones(n,4);

for i =1:n
[X_sol(i,:),fval] = fmincon(@(X)Fig3_solve(X,r(i),lambda_b,lambda_g,gamma,theta_b,theta_g,p,q,theta_bb,theta_gg,r_hat,U),x0,A,b,Aeq,beq,lb,ub,[],options);
end

x_b = X_sol(:,1);
y_b = X_sol(:,2);
x_g = X_sol(:,3);
y_g = X_sol(:,4);


plot(r,x_b,'-',r,y_b,'--',r,x_g,'g--x',r,y_g,'--o','Linewidth',1.5)
legend("\bf x_B","\bf y_B","\bf x_G","\bf y_G",'Location','northwest')
xlabel('sex ratio (ratio of girls to boys), r') 

